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ABSTRACT 

In  order  to  demonstrate  a  promising  technique,  the  control  gain  for 
a  nonlinear,  first  order  system  with  limited  control  and  with  random  input, 
required  to  obtain  a  limited  range  of  state  with  given  likelihood  is  evalu¬ 
ated,  The  proposed  technique  is  based  on  a  series  approximation  for  the 
solution  of  the  partial  differential  equation  for  the  joint  probability 
density  function  of  the  state  variable,  i.e.,  the  Fokker-Pl anck  Equation. 

For  the  simple  illustrative  example,  the  exact  stationary  density  function 
is  found  to  give  control  gains  in  good  agreement  with  the  proposed  series 
approximation  technique.  The  exact  and  approximate  variances  of  state  are 
also  found  in  good  agreement,  but  the  fourth  order  moments  begin  to  show 
some  discrepancy  and  the  detailed  density  functions  do  not  show  good  agree¬ 
ment.  However,  for  many  control  problems  and  for  system  performance  analy¬ 
sis,  accurate  prediction  of  means  and  second  order  moments  is  sufficient, 
as  was  found  for  the  illustrative  example. 

The  series  solution  technique  is  further  demonstrated  by  calculating 
the  transient  state  density  function  with  given  initial  condition  for  the 
same  nonlinear,  first  order  system.  The  results  are  discussed  and  elaborated 
in  order  to  indicate  the  significance  of  this  promising  technique. 
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NOMENCLATURE 

Control  gain  in  system  Equation  (2.17) 
i-th  value  of  a  in  Newton-Raphson  iteration 
Function  defined  in  Appendix  A,  Eq.(A-13) 

Constant  of  integration,  Eq.(2.6) 

Expectation  operation  or  ensemble  average 
Error  function  as  defined  in  Reference  A!  (p.33) 

General  vector  system  function 
Stationary  system  function 

Derivatives  of  f(x)  with  respect  to  control  gain  a 
i-th  component  of  f[x(t),t] 
i-th  coefficient  in  series  expansion  of  f(x) 
i-th  coefficient  in  series  expansion  of  f  (x) 

Constant  driving  noise  gain 
General  driving  noise  gain  matrix 
i , j  — th  element  of  matrix  product  GQG^ 

Function  defined  in  Eq.(2.22) 

Der i vati ve  of  h(a) 

Observation  model 

Hermite  polynomial  of  n-th  order  summarized  in  Appendix  A 
Target  set 

Truncation  value  of  series  approximation 
k-th  moment  about  mean 
Stationary  state  density  function 
Density  function  for  initial  state 
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Density  function  for  state  at  time  t 

Conditional  density  function  for  state  at  time  t  ,  given  the 
density  for  the  initial  state  xQ  at  time  tQ 

Conditional  density  function  for  state  at  time  t  ,  given  the 
observation  set 


P  (^lxu=5)  Conditional  density  function  for  observation  y,  given  the 
y  K  K 


state  x^=§ 

Desired  likelihood  level 

n-th  coefficient  in  series  expansion  of  p(§) 

n-th  coefficient  in  series  expansion  of  px (§, t ;xq, tQ) 

Variance  of  w(t)  in  first  order,  stationary  system 
Covariance  matrix  of  w(t) 

Conditional  density  function  for  state  at  time  t  ,  given  the 
initial  state  is  a  at  time  0 

Saturation  function  defined  in  Eq.(2.l8) 

Matrix  of  spectral  densities  for  state  vector  x(t) 

T  i  me 

Initial  time 

Time  of  k-th  observation 
Time  to  go 

Unit  step  function  defined  on  page  9 
Random  error  in  k-th  observation 
Driving  noise  process 
State  vector 

Conditional  mean  of  x(t),  given 
Initial  state  vector 

Desired  bound  for  first  order  state  variable 
Value  of  state  vector  at  k-th  observation  time 
k-th  observation 
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Set  of  observed  values  up  to  time  t^ 
Realization  variable  for  initial  state 
Launching  window  width 
Kronecker  Delta  function 
Mean  value  of  state  vector 
k-th  moment  about  5=0 
Realization  variable  for  state 
i-th  component  of  § 
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INTRODUCTION 

A  class  of  problems  which  frequently  arises  in  control  engineering  is 
that  of  designing  or  analyzing  a  system  in  which  some  components  are  fixed 
before  the  control  problem  is  attacked,  either  by  the  process  it  is  desired 
to  control  or  by  restrictions  on  space,  weight  and  cost.  It  is  invariably 
found  that  some  of  the  fixed  attributes  of  the  system  are  nonlinear  since  no 
element  of  a  physical  system  can  be  linear  over  an  infinite  range  of  outputs. 

In  the  case  of  marine  vehicle  dynamics,  some  sources  of  nonlinear  effects 
are:  mechanical  control  limits,  hydrodynamic  pitch-yaw  coupling  due  to 
vertical  asymmetry,  temporal  changes  in  wetted  or  submerged  geometry,  hydro- 
dynamic  influences  due  to  proximity  of  the  f ree-surface.  These  attributes 
are  inherently  and  significantly  involved  in  the  dynamics  of  marine  vehicles 
so  that  ignoring  them  by  treating  a  linearized  system  for  control  design  will 
lead  to  vehicles  which  are  ineffectual  or  unstable,  or  exhibit  poor  perform¬ 
ance  qual i ty. 

In  practical  appl ications,  it  is  common  to  find  that  the  input  quantity 
to  a  control  system  is  not  known  exactly,  either  because  it  is  the  result  of 
random  disturbances  acting  on  the  physical  system  or  because  the  input  signal 
results  from  an  imperfect  measurement  of  a  physical  quantity.  Another  source 
of  uncertainty  in  control  analysis  and  design  problems  is  imperfect  knowledge 
of  system  parameters  in  the  final  product.  Recognition  of  these  sources  of 
uncertainty  has  led  to  the  development  of  design  techniques  for  linear  sys¬ 
tems  with  random  inputs  or  uncertain  parameters.  However,  for  nonlinear 
systems  with  random  inputs  or  uncertain  parameters,  the  available  design 
techniques  are  either  cumbersome  and  financially  impractical  or  they  involve 
approximations  which  have  not  been  found  applicable  to  marine  vehicle  dynamics. 
For  example,  the  Monte  Carlo  simulation  technique  requires  the  computation  of 
an  output  sample  from  which  the  desired  operational  statistics  can  be  esti¬ 
mated  for  a  given  design.  Then  control  system  parameters  must  be  adjusted 
and  another  simulation  carried  out.  This  trial-and-error  process  must  then 
be  repeated  until  an  adequate  design  is  found.  Consequently,  for  a  compli¬ 
cated  system  such  as  a  dynamic  marine  vehicle  with  si x-degrees-of-f reedom,  the 
Monte  Carlo  technique  is  economically  not  attractive  as  a  control  design  tool. 
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Equivalent  linearization  techniques  yield  results  which  are  intuitively 
significant  and  reliably  predict  trends  in  performance,  but  are  not  suf¬ 
ficiently  accurate  for  quantitative  control  system  design.  Stochastic 
approximation  techniques  rely  on  the  assumption  that  the  system  response 
time  is  much  larger  or  much  smaller  than  the  range  of  significant  periods  in 
the  input  process,  neither  of  which  is  true  for  many  marine  vehicles. 

Consequently,  in  the  design  or  analysis  of  marine  vehicle  control  sys¬ 
tems,  the  designer  is  faced  with  a  nonlinear  system  with  random  inputs  but 
the  available  tools  are  limited  to  those  for  linear  systems.  This  study  pre¬ 
sents  a  demonstration  of  one  possible  technique  for  stochastic  analysis  of 
nonlinear  systems.  The  demonstration  is  carried  out  for  the  simplest  of  non¬ 
linear  systems  since  the  purpose  here  is  only  to  illustrate  results  and 
feasibility.  Subsequent  studies  must  deal  with  the  practical  problems  of 
implementing  the  technique  and  defining  its  range  of  applicability.  It  should 
be  noted  that  the  rigorous  mathematical  foundation  for  this  technique  is  not 
presently  available.  Therefore  the  mathematical  conditions  for  which  the  pro¬ 
posed  series  solution  converges  and  those  for  which  a  temporally  stable 
approximation  exists  are  not  known.  This  of  course  is  a  disadvantage  but  it 
should  be  remembered  that  engineers  have  used  mathematical  tools  previously 
without  rigorous  mathematical  foundation.  For  example,  Dirac  Delta  functions 
were  used  extensively  before  the  theory  of  distributions  or  generalized  func¬ 
tions  was  developed  by  the  mathematician.  Furthermore,  the  rigorous  foundation 
for  many  of  the  approximate  techniques  which  have  been  attempted  is  just  as 
much  undeveloped  as  that  for  the  technique  proposed  here. 

The  proposed  technique  will  be  demonstrated  for  a  continuous  dynamical 
system  with  an  n-component  state  vector  x(t)  governed  by  a  differential 
system  which  can  be  written  formally  as 

x(t)  =  f[x(t),t]  + G[x(t) ,t]w(t)  ;  t  ^  tQ 

with  initial  condition 

*(t0)  =  *0  ' 

The  driving  noise  w(t)  must  be  a  white,  Gaussian,  zero-mean  vector  with  m- 
components.  The  initial  condition  xq  may  be  a  random  vector  with  density 
P0(§)  not  necessarily  Gaussian  and  must  be  independent  of  w(t)  for  t^  tQ  . 

The  system  function  vector  f[x(t),t]  and  the  driving  noise  gain  matrix 
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G[x(t),t]  must  satisfy  the  continuity  and  other  conditions  stated  by  Jazwinski 
in  his  Theorem  4.5.  Then,  the  conditional  state  density  function  given  the 
density  of  xQ  ,  satisfies 


37  px^,t;xo*  to) 


_  v  o 

"  _i=i 


(Pxfi)  +  7  S 


n 

L 


Lj  t-i  "XT 

=1  j=i  °s 


with  initial  condition 

Px(5>t0;x0,t0)  =  P0<S) 

where  §  is  an  n-component  vector,  fj  denotes  the  i-th  component  of  f[x(t),tj, 
Q  is  the  covariance  matrix  of  w(t)  ,  and  ( GQ.g"^)  - j  denotes  the  i,j-th  element 
of  the  matrix  product.  The  resulting  solution  with  p  ( §) =6( §-xQ)  is  the  funda¬ 
mental  statistical  descriptor  for  the  given  system  and  the  control  analysis 
problem  can  be  solved  once  this  function  is  evaluated. 

In  this  study,  the  state  vector  was  limited  to  one  component  and  the 
system  treated  is  therefore  a  first  order  system.  The  system  function  is 
assumed  t ime- i nvar i ant ,  and  the  driving  noise  gain  is  assumed  constant.  The 
system  equation  then  reduces  to 

x(t)  =  f[x(t)]  +  g  w(t)  ;  t  s  tQ 


with  initial  condition 


x(to)  =  xo 


and  the  partial  differential  equation  for  the  conditional  state  density  re¬ 
duces  to 


37  Px^.tJXo’to)  "  -  If  (Pxf)  +  7  a 


3  PX 


with  initial  condition 

PxU.t0;x0,t0)  =p0(I) 


and  where  the  driving  noise  has  also  been  assumed  stationary.  Under  these 
conditions  and  if  f [ x ( t) J  represents  an  asymptotically  stable  system,  then  a 
stationary  density  p(5)is  assumed  to  exist  for  t  sufficiently  greater  than  tQ. 
This  stationary  density  function  must  satisfy 


1  s 

2  49 


=  0 
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These  equations  for  the  conditional  state  density  Px (§, t ;xQ, tQ)  and  the 
stationary  state  density  p(§)  are  solved  in  this  report  by  means  of  the 
proposed  series  approximation  technique. 

In  the  next  section  of  this  report,  a  very  specific  control  problem 
is  stated  and  solved  using  the  stationary  density.  One  of  the  reasons  for 
treating  the  stationary  first  order  system  is  that  the  above  differential 
equation  for  p(§)  can  be  solved  exactly,  thus  providing  a  basis  for  compari¬ 
son  with  the  approximate  series  solution.  Results  from  the  approximate 
series  solution  are  illustrated  and  discussed.  Then  the  series  approximation 
technique  for  the  transient  conditional  density  is  demonstrated.  The  results 
are  elaborated  and  discussed. 
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THEORETICAL  ANALYSIS 
The  Stationary  Density 

The  system  considered  here  has  scalar  state  variable  x(t)  which 
satisfies  a  system  equation  of  the  form 

x(t)  =  f[x(t)]  +  gw(t)  (2.1) 

where  f(x)  is  the  system  function,  g  is  the  driving  noise  gain  assumed  con¬ 
stant,  and  w(t)  is  a  stationary  white  noise  process  with  normal  probability 
density  with  zero  mean  and  unit  variance.  To  insure  the  existence  of  the 
probability  density  for  the  state  variable  x(t)  and  to  insure  that 
it  approaches  a  stationary  form  when  initial  transients  are  ignorable,  it  is 
assumed  that  f(x)  and  g  are  time-invariant  as  indicated,  that  w(t)  is 

stationary  and  that  there  exist  equilibrium  states  which  are  asymptotically 

2  1 

stable.  These  mathematical  concepts  are  discussed  by  Leipholz  and  Jazwinski 
and  will  not  be  developed  further  here. 

The  stationary  probability  density  function  p(|)  for  the  state  vari¬ 
able  x(t)  must  satisfy  the  stationary  form  of  the  Fokker-Planck  Equation* 
given  by 

7  g2  tf(S)p(§)3  =  o  (2.2) 

as  well  as  the  conditions 


1  im  p(5)  =0 

|5|-°° 

(2.3) 

CO 

J  p(C)d?  =  1  . 

(2.4) 

•00 


In  this  section,  the  exact  solution  of  this  boundary  value  problem  will  be 
derived,  then  a  series  solution  technique  will  be  derived.  These  solutions 
will  be  illustrated  by  application  to  a  specific  system. 

Equation  (2.2)  can  be  integrated  once  to  obtain 

I  g2  cjf  "  f(5)p  =  C1 


5 


R-2033 


where  Cj  is  an  arbitrary  constant  of  integration.  This  equation  is  a  first 
order  ordinary  differential  equation  with  general  solution  given  by 

2C 

P(S)  =  exp[~  J,f(§)d§}|— 4  Jexp{-  J*f  (§) d§}d§  +  C2J  .  (2.5) 

g  9  g 

No  systems  have  yet  been  found  by  this  author  which  admit  asymptotically 
stable  equilibrium  states  and  for  which  the  Cj-term  in  Eq . (2.5)  satisfies 
Eq.(2.3).  Thus,  it  is  presumed  here  that  p(§)  must  be  of  the  form 

P(S)  -  C  exp{”£  Jf  { 5)  d§|  (2.6) 

where  the  arbitrary  constant  C  is  evaluated  from  Eq.(2.4)  as 

CO 

c'1  =  Jexpj^  Jf(5)d§}d§  .  (2.7) 

-00  9 

Equations  (2.6)  and  (2.7)  define  the  exact  solution  for  the  stationary 
probability  density  function  of  the  state  of  the  system  given  by  Eq.(2.1). 
This  general  result  was  obtained  easily  for  the  first  order  system, but  for 
higher  order  systems  the  exact  solution  of  the  corresponding  Fokker-Planck 
Equation  is  not  generally  available.  Consequently,  an  approximate  solution 
technique  is  attempted  here  in  the  hope  that  it  may  be  more  generally  appli¬ 
cable  to  higher  order  systems. 

It  is  assumed  that  the  stationary  probability  density  p(£)  can  be  ex¬ 
panded  in  a  series  of  Hermite  polynomials  Hen(§)  of  order  n,  in  the  form 

P(5)  -  e~§2/22  P  He  (?)  (2.8) 

n=0  n  n 

where  pn;  n=0,l,2,  ...  are  unknown  constants.  Properties  of  the  complete 
orthogonal  set  of  functions  Hen(5)  are  briefly  described  in  Appendix  A. 
Furthermore,  it  is  shown  in  Appendix  B  that 


since  p(§)  must  satisfy  Eq.(2.4).  The  exponential  term  in  Eq.(2.8)  insures 
that  p(§)  will  satisfy  Eq . (2.3)  and  it  remains  to  determine  the  constants 
Pn  such  that  p(§)  satisfies  Eq - ( 2 . 2) . 

Substituting  Eq.(2.8)  into  Eq.(2.2),  multiplying  by  Hem(§)d§  and 
integrating  over  all  5  yields 
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„(5)>5 


-  1  Hem<5>  jr[e"!S/2He„(5)f(i)]'15}  -  0  .  (2.10) 


The  first  integral  can  be  integrated  by  parts  twice,  making  use  of  Eq.(A-9) 
and  (A-l) ,  to  obtai n 


!"•„(»  -4  Le'?a/2He„<s)id5  ■ 0  :  m  - 


®  m!  6„  m^2 

n  ,m-Z 


(2.M) 


The  second  integral  in  Eq.(2.10)  can  be  integrated  by  parts  in  the  same  fashion 
to  obtain 


K<5>-^Le'5a/2lten<5)f(5)]d5-0  J 

-co  dS 


m  =  0 


jHe^^Oe"5  /2Hen(5)f(§)d§;  m^l  . 


It  is  now  assumed  that  the  system  function  can  be  expressed  as  a  Hermite  series 
in  the  form  (A-5) 


f(5)  -  ^  f.He.tS) 

1=0 


(2.12) 


where 


oo  2 

:imr= -  I  f(5)e"§  /z  Hex(S)dg  . 

V2TT  V.  -®  * 


(2.13) 


Equation  (2.12)  is  then  substituted  in  the  above  integral  and  Eq.(A-13)  is  used 


to  obtain 


oo  2 

jHejS)  4|_e_§  /2Hen(§)f(§)]d§  =  0  ;  m=0 


=  -m  V2tt  2j  f«A(A,m-l  ,n)  ;  m&I  .  (2.14) 

t= 0 


With  Eqs.(2.9),  (2.11)  and  (2.14),  Eq.(2.10)  now  gives 
m=  1  : 

nil  P"L  £  ±  J0^A(X,0.0) 


I 
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m  ^  2 ; 


CO  OO-i 

E  p  I  g2mj 6  ,  +  m^  f„A(£,m-l ,n) 

n=1  wnL2  n,m-2  i=Q  l  J 


-Ljig2m!6ft  ,  +  m  D  f  A(£,m-1 ,0)1 

^L2  0,m-2  £=q  * 


But, 

A(X,0,n)=n!6je>n  ; 

A(M,0) 

b£,0  1 

;  A(£,m-1 . 

,0)  «  (m“1)!6x,m- 

1 

so  that 

these 

results  become 

m  =  1  j 

00 

>) 

p  n! f  - - — 

f 

(2.15) 

n=1 

n  n  72W 

0 

and 

m  ^2 ; 

CO 

2 

n=l 

pnLl  ^ ^n,m-2  ' 

00 

je=o  1 

A(j2,m-l  ,n)  1 
(m-1) !  J 

1 

1  2a 
,2  9  6m,2 

+  ^m-lJ 

(2.16) 

If  the 

series 

expansion  for 

p(§)  converges, 

then  i  t 

can  be  truncated 

at 

some  order  M  and  Eqs.(2.15)  and  (2.16)  then  give  a  system  of  M  linear 
non-homogeneous  equations  in  M  unknown  constants  Pn  I  n=1,2,  . ..»  M,  so 
that  the  series  expansion  coefficients  can  be  determined. 

The  exact  and  series  solution  techniques  will  be  demonstrated  by  appli¬ 
cation  to  a  nonlinear  system  with  limited  control.  In  particular,  consider 
the  system  defined  by 

x ( t)  =  X  -  x3  +  a  sat(x)  +  w(t)  (2.17) 

where 

-1  ;  x  <  -l  , 

sat(x)  =  x  ;  I x 1  £  1  ,  (2.18) 

1  ;  x  >  1 

This  system  with  no  control  (a=0)  has  two  asymptotically  stable  equilibrium 
states  at  x  =±l  ,  while  x=0  is  an  unstable  equilibrium  state.  The  control 
function  has  gain  a  for  Ixl^l  and  is  limited  to  ±a  for  1x|  >1  .  The 


8 


R-2033 


objective  of  the  control  in  this  example  is  to  obtain  a  stationary  proba¬ 
bility  pj  that  the  state  x(t)  at  any  time  t  is  between  -x^  and  +x^ 
for  a  given  x^  >0  ,  i.e.,  find  a  such  that 


A1 

1 

■x. 


p(§)d§  =  p, 


(2.19) 


No  attempt  will  be  made  to  argue  that  this  problem  is  related  to  any  specific 
practical  situation  although  it  is  hoped  that  the  interested  reader  will  see 
the  many  practical  aspects  embodied  in  this  simplified  example.  The  purpose 
of  this  application  is  simply  to  demonstrate  the  series  solution  technique 
in  comparison  with  an  exact  solution  of  the  same  problem,  as  well  as  to  in¬ 
dicate  the  usefulness  of  the  probabilistic  design  approach. 


The  system  function  for  this  example  is  seen  to  be 

f (x)  =  x  -  x3+  a  sat(x)  ,  with  g  .  (2.20) 

Substituting  this  function  in  the  general  solution  for  p(§)  given  by  Eq.(2.6) 
leads  to 

P(S)  =  C  exp{  |1  -  |1  +  £|1  u(l-§2)+  a(|§|-  i)u (f-1)} 
with  (2.21) 

C_1  =  2  Jexp{  ^  u(l-f)+  a(  1 1|  -  i)u(§3-1)}  d§ 

where  U(x)  is  the  unit  step  function 


In  the  expression  for  C  '  ,  use  was  made  of  the  symmetry  of  the  integral  of 
an  antisymmetric  system  function. 


In  order  to  find  the  value  of  a  such  that  p(§)  in  Eq.(2.21)  satisfies 
Eq.(2.19)f  a  numerical  iteration  procedure,  the  Newton-Raphson  technique,  was 
used  which  requires  that  p(£)  be  considered  a  function  of  the  parameter  a  , 
so  that  for  the  i-th  estimate  a.  of  the  desired  value,  Eq.(2.19)  becomes 
x, 


1 


h(a|)  =  J  p(§»ai)d§  -  p,  . 


(2.22) 
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Then  the  next  estimate  a.+^  is  obtained  from  the  Taylor  series  expansion 

h(ai+1)  =  h(a.)  +  (a.+1  -  a.)h'(a.)  =  0 
or 

a.+)  -  a.  -  h(a.)/h'(a.)  (2.23) 

where  h'(a)  is  the  derivative  of  h(a)  with  respect  to  a  .  From  Eq - (2.21) 
and  Eq.(2.22),  h'(a)  is  found  to  be 


h'(a)=  J*  P(S){[_4  U(l-§2)+(l|l-  i)U(§3-l)J  +  dS 

-x. 


whe  re 


(2.24) 


I  |£  =  -2C  Jexp{  f-  -  U(l-f)+  a(  |  §1  -  i)U(g3-l)} 

'  “  0 

•  [f- u(i-§2)  +  (151-  i)u(f-i)]  d§  . 


The  solution  to  this  same  problem  is  obtained  via  the  series  solution 
technique  by  substituting  the  system  function  Eq.(2.2D)  into  Eq.(2.13)  and 
using  the  resulting  coefficients  in  Eq .(2.15)  and  Eq.(2.l6)  to  find  the  un¬ 
known  coefficients  pn  for  the  expansion  of  the  density  function  shown  in 
Eq.(2.8).  The  coefficients  for  the  expansion  of  the  system  function  are 
obtained  from  Tables  A-l  and  A-3  as 


f j  =  a  erf 


(_L)  -2 

vr' 


f3  -  -a 


=  -a 


2U1 


-1/2 

-£—  He.  ( 1)  -  1 

3V2 rT 

2e"1/2  He2l-1 

*/2tT  ^X+i)« 


&  =  2,3,  ... 


f2 r  0  ;  *  “  ••• 

Due  to  the  rapid  increase  in  numerical  difficulties  with  increasing  number  of 
terms  in  series  expansions,  it  is  important  to  take  account  of  any  simplifica¬ 
tions  that  can  be  obtained  from  further  analytical  considerations  before  going 
to  Eq.(2.15)  and  Eq.(2.l6),  such  as,  accounting  for  the  symmetries  in  the 
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particular  system  function  being  considered.  In  this  application,  it  ?s 
seen  that  the  system  function  in  Eq.(2.20)  is  antisymmetric  about  x=0,  i.e., 

f(-x)  =  -f(x) 


leading  to  zero  values  for  the  even  order  coefficients  shown  above.  But  if 
f(§)  in  Eq.(2.2)  is  antisymmetric  then  p(|)  must  be  symmetric  about  |=0 
so  that  all  odd  order  coefficients  P2n+j  must  vanish.  These  symmetries  can 
be  accounted  for  by  relabeling  coefficients  in  Eqs.(2.8)and  (2.12)  and  (2.13) 
in  the  form 


p(C)  =  e'?  /2  2J  p  He ,  (§)  ;  p  «  - -  J*  p(§)He  (§)d§ 

n~0  n  2n  n  VH(2n)!-i  2n 


(2.25) 


and 


f(§)  -S  Me  ,(§);  f,  * 1 -  f  f  (§)  e^^He,-..  (?)  dg 

£21+1  £  (2£+l)  !  -t  2£+1 


jfc=0 


(2.26) 


As  before,  Eq.(2.4)  still  requires  pQ  =  l/vTrr  .  Equations  (2.15)  and  (2.16) 
become 


.  v  r  A(2£+l.  2m- 1 ,  2n) 'l 
+  t  l  (2m- 1)!  J 


1  I  1  a  -  ,  jn-1 

Vi 


m  =  1,2, 


(2.27) 


which  again  gives  a  linear  system  of  non-homogeneous  equations  for  the  unknown 
coefficients  pn  ,  once  the  n  and  £  summations  are  truncated. 

In  order  to  apply  Eq.(2.23)  to  the  given  problem  by  means  of  the  series 
solution  technique,  a  representation  for  dp/3a  must  be  developed  so  that 
h,(aj)  can  be  calculated  from 

X1  a 

h'(a)  =  J  |Ip(5)dg  (2.28) 

"X1 


is  considered  a  runct 


where  as  before,  p(g)  . _  _  _ _ ... 

then  Eq.(2.2)  can  be  differentiated  with  respect  to 

_ _  .  _-t-  _  _  I  +■  I  W  * 

2 


a  to  give 


I  g3  |f  -  If  [f(S)  If] =  If  |_fa(s)p(§)_ 


But 


(2.29) 
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where 

f,<»  -  f<{)  • 

In  Eq.(2.29),  the  series  expansion  of  f(§)  is  given  by  Eq.(2.26)  and  that 
for  p(5)  is  given  by  Eq.(2.25)  with  coefficients  obtained  from  Eq.(2.27) 
for  a  given  value  of  a  .  Furthermore,  the  coefficients  f^  are  seen  to  be 
linear  in  the  parameter  a  so  that  the  corresponding  coefficients  in  the 
expansion  of  fo(§)  are  easily  found  to  be 

I  erf  ( 1///2)  ;  £=0 

-He,  (1)/3v5tt5-  ;  im  1  (2.30) 

-2He2 ,  ( 1 )  A/5ne  (2 M  )  J ;  X  2  2 

With  f(§),  fg(l)  and  p(|)  all  known  in  series  form,  dp/da  can  also  be  expressed 
in  series  form  as 


=  e-S2/2  2 


1 


,  PnHe2n(§);  Pn  =  ‘r- 
n=1  </2n  (2n)  !  -» 


p(9e-5V2He2n(!)d5 


(2.31) 


and  the  unknown  coefficients  p^  can  be  found  from  the  differential  Eq.(2.29). 
The  equations  for  p^  become 


=  -2 


Um-l)! 


(2.32) 


which  can  be  solved  numerically  once  the  n  and  l  summations  are  truncated. 
Thus,  the  Newton-Raphson  iteration  technique  in  Eq.(2.23)  can  be  applied  when 
P(§)  is  expressed  in  series  form,  to  determine  the  desired  value  of  the  con¬ 
trol  gain  a  . 
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The  Transient  Density 

The  transient  density  function  Px( §, t ;xq, t  )  can  also  be  evaluated 
by  means  of  the  series  approximation  by  allowing  the  coefficients  to  be  time 
dependent,  i.e., 

2  co 

Px(§.t;xo>t0)  =  e"^  fl  Pn^t;x0’  VHen^  * 


where 


(3.1) 


Pn^^o’V  =1F~7  J  Px^.t;x0,t0)Hen(§)d§. 

a/2it  n:  -oo 


The  conditional  arguments  x  ,t  will  be  dropped  for  clarity.  Equation(2.4) 
still  requires 


pQ(t)  =  ia/K-  . 

Substituting  Eq.(3.1)  into  the  Fokker-Planck  Equation  given  by 

3  .  1  a  d2  d  , 

3t  Px(5>t)  =  2  99  ^  px  ~  Sf  (pxf) 

then  multiplying  by  He  (§)d§  and  integrating  over  all  §  leads  to 


(3.2) 


(3.3) 


S  p„(t)  jV^2  Hen,(l)Hen(5)d{ 

■  -co 


■  Jo  p"(t)  qq\K<5>  I? 

-J  He„(5)  [e-^/2Hen(S)f(S)]d€>  • 

-00 

Using  Eqs.(3.2)  and  (2.12)  and  the  results  in  Appendix  A  gives 
P0(t)  =  0 

CO  00 

Pl(t)  =  n5,  pn(t)n!je50  f*62£+1,n 

(>  (t)  2  \j  qg26  ?  D  f}  M2.**1 ;  m=2 , 3 , . .  • 

m'  '  n=j  L2  ^  n,m-2  £_g  x,  (m-l).1  J 


which  is  a  linear  system  of  first  order  differential  equations.  The  time 
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history  of  the  unknown  coefficients  can  be  evaluated  by  numerical  integra¬ 
tion,  once  the  n  and  SL  summations  are  truncated  and  once  an  initial 
condition  is  specified.  In  general,  the  initial  condition  is  obtained  in 
the  form 

po(0)  =  1/  ,  (3.5a) 

Pn ( 0)  =  — ]z -  J  P  (S)He  (|)d§  ;  n=l  ,2, . . .  (3-5b) 

V2ti  n!  -® 


An  illustrative  calculation  of  the  transient  density  function  is  shown 
in  Figure  5  where  the  initial  condition  is 


Pc(5) 


_  e 


-fnoi 


Ji 


TJOr 


(3.6) 


The  practical  value  of  this  initial  condition  will  be  demonstrated  in  the 
following  section.  Discussion  of  Results.  In  Figure  5,  it  is  seen  that  the 
transient  density  evolves  from  a  fairly  peaked  initial  condition  into  the 
flat  stationary  density  function  shown  in  Figure  3a  for  the  same  value  of  con¬ 
trol  gain. 
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DISCUSSION  OF  RESULTS 
The  Stationary  Density 

The  control  gain  a  required  to  obtain  a  limited  range  of  state  vari¬ 
able  -Xj  <  x(t)  <  Xj  with  0.90  likelihood  was  calculated  for  the  first 
order  system  defined  by  Eq .(2.17)  using  the  exact  stationary  density  function, 
Eq.(2.21)  as  well  as  the  approximate  series  expansion  Eq.  (2.25)  with  coef¬ 
ficients  given  by  the  solution  of  Eq . (2 . 27) -  The  results  obtained  are  compared 
in  Figures  l,  2,  3  and  h. 

The  exact  system  function  given  by  Eq.(2.20)  is  exhibited  in  Figure  1 
for  the  case  a=-3  ,  together  with  truncated  series  approximations  with  5  odd 
terms  (up  to  9th  order  Hermite  polynomial)  and  13  odd  terms  (25th  order  poly¬ 
nomial).  The  coefficients  used  in  the  series  expansion  are  shown  in  Eq.(2.26). 
It  is  seen  in  the  Figure  that  the  exact  and  approximate  system  function  agree 
very  well  for  the  range  of  x  shown  except  near  the  discontinui ty  in  f(x) 
at  x=l  .  At  this  point,  the  polynomial  approximation  is  continuous  and  there¬ 
fore  tends  to  smooth  out  the  discontinuity.  In  subsequent  calculations,  a  nine 
odd  term  series  was  used  (17th  order  polynomial),  which  lies  between  the  two 
approximations  shown  in  Figure  1. 

The  value  of  control  gain  a  required  to  limit  the  state  x(t)  to  the 
range  -x^  <  x(t)  <  Xj  with  0.90  likelihood  is  plotted  in  Figure  2  as  calcu¬ 
lated  by  means  of  the  exact  density  and  by  means  of  the  approximate  series 
expansion  of  the  density.  It  is  seen  that  the  two  solutions  for  the  stated 
control  problem  agree  well.  The  estimated  control  gains  based  on  the  approxi¬ 
mate  density  show  some  scatter  or  irregularity  which  could  be  reduced  by 
improving  the  iteration.  Consequently,  the  proposed  series  approximation 
technique  for  nonlinear  systems  analysis  is  considered  worthy  of  further  devel¬ 
opment  and  study. 

It  should  be  noted  here  that  the  level  of  effort  required  to  obtain  the 
results  shown  in  Figure  I  by  means  of  the  Monte  Carlo  simulation  technique 
(the  best  known  technique  when  the  exact  solution  is  not  obtainable)  is  orders 
of  magnitude  greater  than  that  utilized  here  in  applying  the  series  approxi¬ 
mation  technique  to  derive  the  density  function.  The  required  steps  would  be 
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as  follows.  A  simulator  or  model  of  the  system  defined  by  Eq .(2.17)  would  be 
developed  and  subjected  to  a  sample  input  process  w(t)  ,  for  a  given  control 
gain  a  .  The  output  sample  would  be  analyzed  statistically  in  order  to 
estimate  the  probability  of  the  event  {-x^ <  x(t) <  x^ \  as  well  as  the  confidence 
of  this  estimate.  This  process  would  then  be  repeated  for  a  new  value  of  a 
until  the  desired  1 ikel ihood  was  obtained,  with  appropriate  confidence.  The 
only  way  to  increase  confidence  in  the  estimated  likelihood  is  to  increase  the 
sample  length  of  the  simulation  or  model  test.  Finally,  all  the  above  steps 
must  be  repeated  for  each  Xj  shown  in  Figure  1.  On  the  other  hand,  the 
series  approximation  technique  requires  one  solution  of  a  linear  algebraic 
system  of  equations  for  each  iteration  step  and  the  results  shown  in  Figure  2 
were  obtained  with  2  to  3  iterations. 

Some  calculated  density  functions  are  shown  in  Figure  3a,  b  and  c,  using 
the  exact  and  approximate  solutions,  for  various  values  of  x^  .  When  the  de¬ 
sired  range  x^  is  large  (x^=1.5  and  1 .35 , F i g . 3a) ,  the  exp(-§4)  behavior 
for  5  large  in  the  exact  solution  Eq . (2.21)  produces  a  rapid  fall-off  of 
p(§)  which  the  series  solution  attempts  to  follow.  However,  the  approximate 
solution  behaves  like  exp(-§2)  times  a  polynomial  so  that  in  attempting  to 
produce  the  more  rapid  fall-off,  poor  agreement  for  small  §  results.  On  the 
other  hand,  for  small  values  of  X|  (0.60  and  0.75  in  Fig. 3c),  the  large  peak 
at  §=0  leads  to  poor  agreement  between  exact  and  approximate  solutions. 

In  between  these  extremes,  some  good  approximations  for  the  exact  density 
function  are  shown.  In  particular,  see  the  result  for  x^=0.90  in  Figure  3c. 
It  should  be  kept  in  mind  here  that  good  estimates  of  the  desired  control 
gain  a  were  obtained  for  all  cases  shown  in  Figure  3  so  that  a  solution 
technique  which  duplicates  the  exact  density  function  in  every  detail  is  not 
always  required.  Methods  for  improving  the  convergence  of  the  series  solution 
technique  and  for  handling  the  numerical  difficulties  of  large  order  systems 
must  be  developed  during  the  subsequent  evolution  of  the  proposed  technique. 

A  further  comparison  between  the  exact  and  approximate  solution  is 
illustrated  in  Table  1  which  shows  corresponding  values  of  the  second  order 
moment  M2=ox2  ,  the  variance,  and  the  ratio  M^/M^  ,  the  kurtosis,  which  is 
a  measure  of  the  flatness  or  peakedness  of  the  density  function.  From  Eq.(B-3) 
It  IS  seen  that  = 3  for  a  normal  density.  The  comparison  between  exact 

and  approximate  variance  is  seen  to  be  very  good  for  all  cases  except  XpO.60 
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where  the  exact  solution  becomes  very  peaked  at  §=0  .  The  discrepancy  in 
density  function  in  most  cases  does  not  effect  the  accuracy  of  the  series 
solution  for  second  order  statistics  just  as  for  estimating  the  control 
gain  a  .  The  comparison  of  kurtosis  is  seen  to  be  very  good  except  for  the 
last  two  values  (xj  =  0.75,  0.60).  Thus,  the  accuracy  of  the  series  solu¬ 
tion  diminishes  for  predicting  higher  order  moments  as  should  be  expected. 

For  performance  analysis  and  for  control  system  design,  estimates  of  the 
mean  and  variance  are  the  most  important  statistics  and  these  are  reliably 
predicted  here. 

The  series  for  PX(S)  was  truncated  at  M=8  in  all  the  results 
shown  in  Figures  2  and  3  and  in  Table  1.  The  effect  of  decreasing  the  num¬ 
ber  of  terms  in  the  series  expansion  is  shown  in  Figure  k.  It  is  seen  that 
the  agreement  between  exact  and  approximate  solutions  diminishes  as  the  series 
is  truncated  at  lower  order.  Solutions  attempted  with  larger  numbers  of  terms 
showed  signs  of  numerical  problems  which  were  beyond  the  scope  of  this  initial 
s  tudy . 


The  Transient  Density 


Utilization  of  the  series  approximation  technique  to  solve  the  Fokker- 
Planck  Equation  has  been  demonstrated  in  Figure  5  for  a  nonlinear  first  order 
system  with  limited  control.  Suppose  that  this  technique  could  yield  a  means  of 
approximating  the  density  function  rx(§,t|«)  defined  as  the  transient 
density  for  the  given  system  with  initial  condition 


rx(S,C|a)  =  6(S-a) 

where  a  is  any  given  initial  state  and  where  a  time  invariant  system  has 

been  assumed  so  that  we  can  taka  t  =0  for  simplicity.  It  should  be  noted 

o 

that 


1  im 


-r/2  ar 


V2Ttc 


=  6(1) 


so  that  the  results  shown  in  Figure  5  are  an  approximation  for  r^ ( § , t | 0) .  In 
any  event,  rx(5,t|c*)  is  the  conditional  transition  density  function  for  the 
system  from  which  any  first  or  second  order  statistic  can  be  calculated. 


The  joint  state  dens' ty, for  example,  is  given  by 
Px( t+T ;cy,  t)  =  ^(S.Tlo)  Px(Q.t)  ,  T  2  0 
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I  * 


from  which  the  power  spectral  density  for  stationary  operating  condition 
Sxx(uu)  can  be  calculated  by  taking  px(o,t)  to  be  the  stationary  density 
Px(a)  and  using 

SxXW  =  JdT  e"llWT  JJd5  da  rx(§,T|a)  px(a) 

-00  -00  ^ 

The  first  order  state  density  is  given  by 

00 

Px(§.t+T)  -J*  p  (§,t+T;a,T)da( 


=  J*  rx(§,rja)  Px(a,t)dof 
-00 

which  is  the  Chapman-Kolmogrov  Equation.  From  the  first  order  state  density, 
the  mean  p,  (t)  and  variance  a  2(t)  can  be  calculated  by  means  of 

X  X 

nx(t)  =  J  5  Px(g,t)d5 

-00 

<(t)  -  J(g-v.x)8  px(g,t)d|  . 

-00 

As  an  illustration  of  how  these  results  can  be  utilized  for  performance 
analysis,  suppose  it  is  desired  that  the  state  of  the  system  after  time  T  be  i n  an 
interval  |^_  .  The  probability  of  obtaining  the  desired  state,  when  the  initial 
state  is  a  random  variable  xq  with  density  PQ(a)  ,  is  given  by 

Probix(T)elTl  =  J  px(?,T)d§ 

!T 

CO 

=  J  J  rxU,T|or)p0(a)da'  . 

IT  - 

The  series  approximation  technique  illustrated  herein  shows  promise  of  pro¬ 
viding  a  relation  for  the  fundamental  statistical  descriptor  of  the  system 
>‘X(§,T|Q')  in  terms  of  system  parameters  such  as  the  control  gain  a  ,  and  in 
terms  of  time  T.  Supposing  further  that  the  initial  density  is  of  the  form 


»(“)  =  { 


1  /2(3  ;  lad  <  0 
0  ;  I  o'  I  >  p 


then  the  designer  could  find  the  combinations  of  launching  window  width  P, 
time  to  go  T  and  control  gain  a  ,  which  would  yield  an  acceptable  risk  of 
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missing  the  target  set  1^  given  by  1-Prob jx(T)  el^.|  . 

Finally,  the  statistical  estimation  of  state  of  a  nonlinear  system 
driven  by  noise,  from  observations  of  state  corrupted  by  noise,  can  be 
carried  out  with  the  series  approximation  technique  outlined  here.  The 
problem  of  parameter  identification  falls  in  this  category.  In  particular, 
suppose  that  the  state  x(t)  of  a  system  is  modeled  by  an  equation  of  the 
form 

x(t)  =  f[x(t),t]  +  G[x(t) ,t]w(t) 

t  ^  t  ;  x(t  )  =  x_ 
o  o'  o 

where  w(t)  is  a  white,  normal,  zero  mean  process,  where  xq  may  be  a  random 
vector  which  is  independent  of  w(t),  and  that  the  observation  model  is  given 
by 

~  ^xk’  +  vk  ’  •  •• 


where  v^  are  independent,  normal,  zero  mean  random  vectors  which  are  in¬ 
dependent  of  xq  and  w(t)  .  The  statistical  problem  is  to  derive  an 
estimate  x(t)  for  the  state  based  on  measurements 

vk  -  bpy2 . yk*  ;  *k  •=  1  s  (k+i  • 


One  good  estimator  is  the  conditional  mean  defined  as 

00 

x(t)  =  E|x(t)|Ykh  J*§  Px(§»t|Y1)d§  . 

•00 

It  has  been  shown^  that  between  observations  the  conditional  density  of  func¬ 
tion  px ( 5 , t I Y^)  must  satisfy  the  Fokker-Planck  Equation 


St  PX  +  2  j?,  L(GQGT)ijPx] 


for  t^  <  t  <  t^+j  with  initial  condition 


Px(§.tklvk) 


At  an  observation  the  conditional  density  function  is  updated  by  means  of 
Bayes's  Theorem  in  the  form 
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px(S<tklV 


Jpy(yk|xk-il)  Px(’l.*kIVi)d11 


which  relates  the  new  initial  condition  to  the  density  at  the  end  of  the 
previous  interval  between  measurements.  Thus,  the  required  conditional 
density,  and  thereby  x(t)  ,  can  be  evaluated  for  all  t  by  alternately 
using  the  Fokker-Planck  Equation  and  Bayes's  Theorem.  The  series  approxi¬ 
mation  provides  the  appropriate  means  for  solving  the  Fokker-Planck  Equation. 


R-2033 


CONCLUSIONS 

For  effective  design  of  nonlinear  systems  subjected  to  random  inputs, 
a  stochastic  analysis  technique  is  needed  in  order  to  estimate  statistics 
of  output  or  performance  for  given  input  statistics.  A  technique  for  pro¬ 
viding  this  relationship  is  illustrated  herein,  based  on  a  series  approxi¬ 
mation  of  the  solution  to  the  governing  equation  for  the  fundamental 
statistical  quantity,  the  transient  conditional  state  density  function. 

Using  the  governing  equation,  the  Fokker-Planck  Equation,  a  stationary 
control  problem  for  a  first  order  nonlinear  system  with  limited  control,  was 
stated  and  solved  exactly  and  by  means  of  the  series  approximation  technique. 
The  exact  and  approximate  control  gain  were  found  to  be  in  good  agreement  as 
were  the  variances  of  response.  The  fourth  order  response  moment  began  to 
show  some  discrepancy  between  the  exact  and  approximate  solutions  and  the 
detailed  density  functions  were  in  poor  agreement,  but  the  stated  problem 
did  not  require  detailed  agreement  as  is  the  case  for  many  design  and 
analysis  problems. 

Also,  based  on  the  governing  equation,  the  series  solution  technique 
for  evaluating  the  transient  conditional  density  function  was  demonstrated 
for  the  same  first  order  nonlinear  system  and  control.  Then,  in  discussing 
these  results,  it  is  briefly  shown  how  the  transient  conditional  density 
function  can  be  used  to  analyze  problems  in  systems  analysis,  control  design 
and  state  estimation. 
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TABLE  1 

COMPARISON  OF  MOMENTS 


*1 

VARIANCE  (M2) 

KURTOSIS  (M4/M22) 

Exact 

Approx. 

Error 

Exact 

Approx. 

Error 

1.50 

0.921 

0.926 

-0.5% 

1.99 

1.97 

1.0% 

1.35 

0.705 

0.725 

-2.8 

2.22 

2.18 

1.8 

1.20 

0.541 

0.522 

3.6 

2.45 

2.50 

-2.0 

1.05 

0.411 

0.402 

2.2 

2.65 

2.67 

-0.7 

0.90 

0.302 

0.300 

0.7 

2.72 

2.84 

-4.2 

0.75 

0.209 

0.207 

1 .0 

2.92 

3.16 

-7.6 

0.60 

0.133 

0.120 

10.8 

2.95 

6.32 

-53.3 
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FIG.  I.  SERIES  APPROXIMATION  FOR  SYSTEM  FUNCTION 
WITH  DISCONTINUOUS  SLOPE 
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90%  INTERVAL  WIDTH 


x. 


2.  LIMITED  CONTROL  CONSTANT  REQUIRED  TO  OBTAIN 
90  PER  CENT  UKLIHOOD  THAT  STATE  IS  BOUNDED  BY  X  , 
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FIG.  3a.  COMPARISON  OF  EXACT  AND  APPROX  IMATE  STATE  DENSITIES 


FIG.  3C.  COMPARISON  OF  EXACT  AND  APPROXIMATE  STATE  DENSITIES 


FIG. 4.  EFFECT  OF  SERIES  TRUNCATION  ON  APPROXIMATE  DENSITY  FOR 
X.  =  1.35 
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APPENDIX  A 

HERMITE  POLYNOMIAL  -  SUMMARY 


The  following  summary  of  properties  of  the  Hermite  polynomial  Hen(x) 

of  degree  n  and  argument  x  ,  follows  the  more  detailed  and  broader  de- 

A1 

velopment  of  Abramowitz  and  Stegun.  This  set  of  functions  satisfies  the 
orthogonality  condition 

fe  X  2He  (x)He  (x)dx  =  */2n  n!  6  (A-l) 

J  n  '  m  m,n  ' 

7 


where  the  order  is  a  non-negative  integer,  the  argument  is  a  real  variable 

and  6  is  the  Kronecker  delta  function  defined  as 
m,n 


r  0  ;  m/n 
^  1  ;  m=n 


(A-2) 


The  Hermite  polynomial  is  further  defined  by  the  finite  series  expansion 


[n/2]  ,  lXm  n-2m 

He  (x)  =  n!  2]  ^ - - 

m=0  m!  2m(n-2m)  J 


(A-3) 


where  the  notation  [n/2]  denotes  the  largest  integer  not  greater  than  n/2. 

A  brief  list  of  coefficients  for  the  Hermite  polynomial  is  shown  in  Table  A-l. 

Using  the  orthogonality  condition,  Eq.(A-l),  a  function  f(x)  for 
which  the  expansion  exists,  can  be  expressed  as  an  infinite  series  in  two 
ways,  namely,  either 

2  .  co  _  co 

f(x)  =  e  X  2j  f  He  (x)  ;  f  =  -  ff(x)He  (x)dx  (A-4) 

"=0  n  n  n  V2W  n!  -i  n 

or 

oo  oo  2 

f(x)  =  2]  f  He  (x)  ;  f  =  - - —  f  e  X  ^2f(x)He  (x)dx  .  (A-5) 

n=0  n  n'^  n  ^  n,  j, 


The  choice  of  which  form  to  use  depends  mainly  on  the  behavior  of  f(x)  as 
Ixl  -*  •  .  If  the  integral  defining  fR  in  Eq.(A-4)  exists,  then  either  form 
may  be  utilized  but  if  this  integral  does  not  exist,  then  Eq.(A-5)  must  be 
used.  The  coefficients  for  the  expansion  of  several  functions  are  shown  in 
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Tables  A-2  and  A-3,  using  Eqs.(A-4)  and  (A-5) ,  respectively. 

From  Eq. (A-3) >  it  is  seen  that 

Hen(-x)  =  (-1)n  Hen(x)  (A-6) 

that  is,  even  order  Hermite  polynomials  are  symmetric  functions  of  x  while 
odd  orders  are  antisymmetric.  The  Hermite  polynomial  of  n**1  order  satisfies 
the  ordinary  differential  equation 

j2  j 

Hen(x)  “  x  dx  HenM  +  nHen(x)  =  0  ,  (A-7) 

it  satisfies  the  recurrence  relation  with  respect  to  n  given  by 

Hen+J(x)  =  xHen(x)  -  nHe^fx)  ,  (A-8) 


its  derivative  with  respect  to  x  is  given  by 


-r-  He  (x)  =  nHe  .  (x) 
dx  n  '  n-1'  ' 


(A-9) 


xz  ™  z3 /2 

and  its  generating  function  is  e  so  that 


3/o  ®  n 

,*2-z  n  .  s  i_  He  (x) 

n=0  n*  n 


(A- 10) 


and 


u  (  \  dn  xz-z3/2 
He  (x)  =  -  e 

n  dz" 


(A-ll) 


z=0 


Rodrigues' 


Formula  for  the  Hermite  polynomial  of  nt*1  order  becomes 


Hen(x)  =  (-l)n 


dn  -x3/2 


(A- 1 2) 


Several  definite  integrals  involving  Hermite  polynomials  are  shown  in 
Table  A-4.  The  notation  in  Table  A-4  follows  that  of  Abramowitz  and 
Stegun.^' 
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Using  the  generating  function  in  the  form  of  Eq.(A-ll),  it  can  be 
shown  that 

eo  -t2 /2 

J  - — ~  He.(t)He.(t)He  (t)dt  =  A(i,j, 

-«  V2TT  '  J  k 

_  i!  ii  k! 

-  ,  (±^ .  . 

=  0  ;  otherwise.  (A—  7 3) 


k) 


if 


i+j+k  is  even  and 
i^j+k  and 
j^k+i  and 
k^i+j  , 
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TABLE  A-l 


LIST  OF  COEFFICIENTS  FOR  Hen(x)  and  x 


Hen(x) 

[n/2] 

=  E  H)m 

m=0 

h(n)xn-2m 
m  * 

-C 

1 — » 

Csl  O 

1 _ 1 

u 

c 

X 

^  He  ,  (x) 
m  n-2mv  ' 

I 

Table  of 

h<"> 

m 

1 

1 

! 

m= 

n 

0 

1 

2 

3 

4 

5 

6 

i 

7  1 

1 

0 

1 

2 

1 

1 

1 

1 

3 

1 

3 

4 

1 

6 

3 

5 

1 

10 

15 

6 

1 

15 

45 

15 

: 

7 

1 

21 

105 

105 

8 

1 

28 

210 

420 

105 

' 

9 

1 

36 

378 

1260 

945 

i 

10 

1 

45 

630 

3150 

4725 

945 

! 

11 

1 

55 

990 

6930 

17325 

10395 

j 

12 

1 

66 

1485 

13860 

51975 

62370 

10395 

13 

1 

78 

2145 

25740 

135135 

270270 

135135 

• 

14 

1 

91 

3003 

45045 

315315 

756756 

945945 

135135 

15 

1 

105 

4095 

75075 

675675 

I891890 

4729725 

2027025 

16 

1 

120 

5460 

120120 

1351350 

5045040 

18918900 

16216200 

17 

1 

136 

7140 

185640 

2552550 

12252240 

64324260 

91891800 

18 

1 

153 

9180 

278460 

4594590 

27567540 

192972780 

413513100 

19 

1 

171  11628 

406980 

7936110 

58198140 

523783260 

1571349780  j 

20 

1 

190  14535 

581400 

13226850 

116396280 

1309458150 

5237832600 

m= 

8 

9 

10 

n 

16 

2027025 

17 

34459425 

18 

310134825 

34459425 

19 

1473140419 

654729075 

20 

7365702095 

6547290750 

654729075 

h<">  -  — 

n! 

M  m!2m(n-2m)! 


i 
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TABLE  A-2 

HERHITE  SERIES  EXPANSIONS 


,,  \  -x2/2 

f(x)  =  e 


U  fnHen(x)  ;  */2n  nif  =  J  f(x)Hen( 
n-0  -os 


f(x) 

1  ;  Ixl  <  1 
0  ;  Ixl  >  I 

e-(x-ta,)2/2cr2 
a/2TT  a 


f 

n 


0  ;  n  odd 

n+T  Hen+1 ( I)  I  n  even 


■  I 

n  -V 


.[n/2] 

m=0 


lill 


.  ~m  2m 

m!2  (i. 


j-n^rnj 

2 

k=0 


NOTES: 

sat(x) 


-1;  x  <  -1 
x;  Ixl  <  1 
+1 ;  x  >  1 


'2  k 


J_  J  §2k  e"^72 

*/2tt 
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TABLE  A-3 

HERMITE  SERIES  EXPANSIONS 

®  eo  2 

f(x)  =  n?0fnHen(*)  ;  ^2TT  n»  f n  -  J*  e"X  /Zf(x)Hen(x)dx 


Ii*L 


erf\~-) 


sat(i) 


2  (-1/3) (n"})/2 

V3n  n(^~ )] 


-2e~bS/2Hen_2(b) 
*/2n  b  n  i 


;  n  even 


n  odd 


;  n  even 


;  n=3,5, 
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TABLE  A-4 


i: 
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Finally,  the  standard  normal  density  is  obtained  by  setting  jx=0  and 
o=l  in  Eq.(B-l),  leading  to 


P(5)  - 


e-S8/2 

72TT 


(B-7) 


which  is  seen  to  be  proportional  to  the  weighting  function  for  the  Hermite 
polynomial  shown  in  Eq.(A-l). 


The  Hermite  series  expansion  of  the  normal  density  is  given  by 

P(S)  -  e_§2/2E if  He  (5)  .  (B-8) 

n=0  n  n 

Multiplying  by  Hem(§)  and  integrating  leads  to 

Jp(S)He  <|)d;  -S  f„  J  (5)d| 

-CO  11  -CO 

so  that  the  orthogonality  condition  in  Eq.(A-l)  gives 


1 


fm  ^  ml 


J  p(5)Hem(S)dg 


Introducing  Eqs.(B-l)  and  (A-3)  gives 

f  .  -l1%2] _ ”r  .-(si*)a/2osdj 

n!2  (m-2n) :  -00 

where  the  integral  is  seen  to  be  the  m-Zn^  moment.  Thus, 


f 


m 


[m/2] 

-L. 

tfn  n=0 


n!2n(m-2n) ! 


(B-9) 


which  defines  the  Hermite  series  coefficients  in  terms  of  the  moments  of 
the  normal  density  function. 


Multiplying  both  sides  of  Eq.(B-8)  by  HeQ(§)=l  and  integrating 

gives 

CO  QO  00  2 

Jp(S)d5=^  f  J  e“5  /2  He  (§)He  (§)d§  . 

-«  n=0  n.»  n  0 

But  the  area  under  any  density  function  is  unity  and  the  orthogonality 


C 
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condition  brings  this  result  to  the  form 

fQ  =  1/VH  (B-10) 

which  agrees  with  Eq.(B-9)  and  demonstrates  that  the  only  term  in  the  Hermite 
series  expansion  which  contributes  to  area  is  the  zeroth-order  term. 

Finally,  combining  Eqs.(B-2),  (B-8)  and  (B-9)  leads  to  the  final  rela¬ 
tionships  between  the  moments  and  the  expansion  coefficients: 


=  \/2tt  f 


1 


=  V2tt 


[m/2] 

!  E 

j2=0 


m-2X‘ 


nln 


( B  —  1  1 ) 


The  central  moments  can  then  be  formed  from  Eq.(B-5)  or  from  the  recursion 
relationship  in  Eqs.(B-3).  A  number  of  explicit  expressions  for  expansion 
coefficients  and  moments  are  given  in  Table  Bl. 


After  several  pages  of  algebraic  manipulation,  it  can  be  shown  that 
the  Hermite  expansion  coefficients  can  be  written  in  terms  of  the  mean  n 
and  variance  a3,  in  the  form 


:  _  J_  >"  (2n) ! 

2"  "VliteO  <**)  !<"-*)! 


^ ) 


n-X 


;  n=0 , I , . . . 


(B-12) 


f2n+l 


-  J Xs  (2n) ! 
y2TT  X=0  (2X+1)  J(n-X)! 


n=0 


With  these  results  it  is  then  found  that  for  the  special  case  p.=0  ,  the 
expansion  of  the  normal  density  becomes 


P(S)  - 


e-§3/2o3 
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Jo  £  (49 
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which  is  seen  to  be  symmetric  in  5  as  required.  Furthermore,  for  the 
special  case  ct=1 ,  Eq.(B-12)  gives 


P(S) 
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while  for  the  case  when  p=0  and  <j=|  ,  the  expansion  reduces  to 

-%zfl  -S2/2 

P(C)  -  *-=7-  =  “H —  He  (?)  .  (B-15) 

V2tt  V2n  ° 

Equations  ( B- 1 3)  and  (B- 1 4)  clearly  demonstrate  that  the  number  of  terms 
required  in  the  expansion  to  obtain  a  desired  accuracy  increases  as  o2 
deviates  further  from  unity  or  as  (i  deviates  further  from  0  . 


The  general  relationships  between  the  probability  density  and  distri¬ 
bution  functions  are 

d  5 

P(5)-||P(S)  ,  P(§)  p(S'W  .  (B-16) 

Thus,  the  normal  distribution  function  for  a  random  variable  with  mean  M* 
and  variance  cr2  is  given  by 

S  «“(5,“**)3/2oa 

-  S  - — ~ - d§‘ 

-00  «/2n  a 


where 


•  H'  *  "fe)] 


erf(z)  =  J*  e-t  dt 

«/TT  o 


(B-17) 


following  Reference  Al.  Integrating  the  expansion  of  the  density  function 
given  in  Eq.(B-8),  term  by  term  yields 


R-2033 


TABLE  B 1 

SUMMARY  OF  EXPLICIT  EQUATIONS  FOR  COEFFICIENTS  AND  MOMENT 
FOR  THE  NORMAL  PROBABILITY  DENSITY  FUNCTION 

^  =  Mean  o2=  Variance 

fn  =  1A/2W 

v  2  co 

f  =  p./n/2tt  p(C)  =  e”5  /2S  fnHen(§) 

•  n=0 

f2  *  (m-2  +  o2  -  1)  /2\/2TT 
f3  -  nip2  +  3 (o2  -  1)J/6s/2tt 
f4  =  [p4  +  -  1)  +  Vf?  -  1)2]/24v/H 

f5  =  io,[p4+  10 ^(o2  -  1)  +  ^(a8-  1)2]/12072tt 

f0 

-  ^TT  fj 

«  Ji*  (2f2  +  fQ) 

P3  =  *JTn  (6f ^  +  3f,) 
p4  =  72W  (24f4  +  12f2  +  3f0) 

P5  =  (120f5  +  60f3  +  15^1) 

1 
0 

(2f2  -  f,  +  f0) 

Jin  (6f3  -  6f]f2  +  2f2) 

</2W  [24f4-24f,f3  +  12f2(f,2+fQ)-  3ff  -  6f ,2  +  3fQ] 
V2Tt[l20f5-  120f ? f4+  60f3(f,2+fQ)-  20f2f|(fJs+3f0)+  20f?3 
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